jacobian <- function(x, prods, T, nn, v, Y, K, x_char, prodsMarket, numProdsTotal, marketStarts, marketEnds, oo, sharesum, marketForProducts, nD) { 

x<-as.matrix(x)

theta1 <- x[1:(K+1), 1]     
theta2 <- x[((K+1)+1):((nD+1)*(K+1)),]    
delta = x[((nD+1)*(K+1)+1):((nD+1)*(K+1)+numProdsTotal), 1]               
g = x[((nD+1)*(K+1)+numProdsTotal+1):nrow(x), 1]  
g <- as.matrix(g) 
  
#theta1 <- x[1:(K+1), 1]                           
#theta2_mat <- matrix(0, K+1, nD)
#for (dd in 1:nD) {
#  theta2_mat[,dd] = x[(dd*(K+1)+1):((dd+1)*(K+1)), 1]
#}
#delta = x[((nD+1)*(K+1)+1):((nD+1)*(K+1)+numProdsTotal), 1]               
#g = x[((nD+1)*(K+1)+numProdsTotal+1):nrow(x), 1]  
#g <- as.matrix(g) 
  
#v_mat <- array(0,dim=c(K+1,ncol(v),nD))
#for (dd in 1:nD) {
#  v_mat[, ,dd] = v[((dd-1)*(K+1)+1):(dd*(K+1)),]
#}


expmu = matrix(0,nrow(x_char),ncol(v))
for (t in 1:T) { 
expmu[marketStarts[t]:marketEnds[t],] = exp(x_char[marketStarts[t]:marketEnds[t],]%*%diag(theta2[1:(K+1)])%*%v[1:(K+1),,t])
for (dd in 2:nD) {
expmu[marketStarts[t]:marketEnds[t],] = expmu[marketStarts[t]:marketEnds[t],] * exp(x_char[marketStarts[t]:marketEnds[t],]%*%diag(theta2[((K+1)+1):(2*(K+1))])%*%v[((dd-1)*(K+1)+1):(dd*(K+1)),,t])
}
}

# v
#expmu<-exp(x_char%*%diag(theta2[1:(K+1)])%*%v[1:(K+1),])

# Income
#for (tt in 1:T) { 
#  expmu[marketStarts[tt]:marketEnds[tt],] = expmu[marketStarts[tt]:marketEnds[tt],] * exp(x_char[marketStarts[tt]:marketEnds[tt],]%*%diag(theta2[((K+1)+1):(2*(K+1))])%*%Y[((tt-1)*(K+1)+1):(tt*(K+1)),])
#}

#expmu <- exp(x_char%*%diag(theta2_mat[,1])%*%v_mat[, , 1])
#for (dd in 2:nD) {    
#  expmu = expmu * exp(x_char%*%diag(theta2_mat[,dd])%*%v_mat[, , dd])
#}

expmeanval = as.matrix(exp(delta))
simEstShare<-ind_shnormMPEC(expmeanval,expmu, oo, sharesum, marketForProducts)
simShare<-simEstShare[,1:(ncol(simEstShare)-1)]
EstShare<-simEstShare[,ncol(simEstShare)]
EstShare<-as.matrix(EstShare)


ooo<-matrix(1,1,K+1)
  
# Evaluate the Gradients
dSdtheta2<-matrix(0,numProdsTotal,(nD*(K+1)))
dSddeltaDIAG<-matrix(0,numProdsTotal,prods)
dSddelta<-matrix(0,numProdsTotal,numProdsTotal)
  

for (t in 1:T) {
  index<-marketStarts[t]:marketEnds[t]
  ooo1<-matrix(1,prodsMarket[t],1)
  index_start <- marketStarts[t]
  index_end <- marketEnds[t]
  
  dSddeltaDIAG_old <- matrix(0, length(index),prodsMarket[t])
  dSddeltaDIAG_temp = dc_1_cpp (simShare, index_start, index_end, dSddeltaDIAG_old, nn )
  dSddeltaDIAG[index,1:prodsMarket[t]] = dSddeltaDIAG_temp
  dSddelta[index,index] = dSddeltaDIAG[index,1:prodsMarket[t]]
  
  dSdtheta2 = dc_2_cpp_new(dSdtheta2, simShare, ooo, ooo1, v[,,t], x_char, nn, index_start, index_end, nD, K, t) 
  #dSdtheta2 = dc_2_cpp(dSdtheta2, simShare, ooo, ooo1, v, x_char, nn, index_start, index_end, nD, K) 
  
}

#for (t in 1:T) {
#  index<-marketStarts[t]:marketEnds[t]
#  ooo1<-matrix(1,prodsMarket[t],1)
  
#  for (rr in 1:nn) {
#          if (length(simShare[index,rr])==1) {
#          dSddeltaDIAG[index,1:prodsMarket[t]] = dSddeltaDIAG[index,1:prodsMarket[t]] + (simShare[index,rr] - simShare[index,rr]*t(simShare[index,rr]))/nn    
#          } else {
#          dSddeltaDIAG[index,1:prodsMarket[t]] = dSddeltaDIAG[index,1:prodsMarket[t]] + (diag(simShare[index,rr]) - simShare[index,rr]%*%t(simShare[index,rr]))/nn  
#          }
#          # v 
#          dSdtheta2[index,1:(K+1)] = dSdtheta2[index, 1:(K+1)] + (simShare[index,rr]%*%ooo)*(ooo1%*%t(v[1:(K+1),rr,t]))*( x_char[index,] - (ooo1%*%(t(simShare[index,rr])%*%x_char[index,])))/nn
#          # Income
#          dSdtheta2[index,((K+1)+1):(2*(K+1))] = dSdtheta2[index, ((K+1)+1):(2*(K+1))] + (simShare[index,rr]%*%ooo)*(ooo1%*%t(v[((K+1)+1):(2*(K+1)),rr,t]))*( x_char[index,] - (ooo1%*%(t(simShare[index,rr])%*%x_char[index,])))/nn
#          
#          #for (dd in 1:nD) {
#          #dSdtheta2[index,((dd-1)*(K+1)+1):(dd*(K+1))] = dSdtheta2[index,((dd-1)*(K+1)+1):(dd*(K+1))] + (simShare[index,rr]%*%ooo)*(ooo1%*%t(v_mat[,rr,dd]))*( x_char[index,] - (ooo1%*%(t(simShare[index,rr])%*%x_char[index,])))/nn
#          #}
#  }
#  dSddelta[index,index] = dSddeltaDIAG[index,1:prodsMarket[t]]
#}

Ddelta = -solve(dSddelta)%*%dSdtheta2


return(Ddelta)
  
}  